using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// The Gauss-Kronrod rule 15 for infinity integrals.
    /// </summary>
    [Serializable]
    public static class Dqk15iClass
    {
        /// <summary>
        /// The Gauss-Kronrod rule 15 for infinity integrals.
        /// </summary>
        /// <param name="f">The function.</param>
        /// <param name="boun">The bound of the integral.</param>
        /// <param name="inf">The value which describes the infinity..</param>
        /// <param name="a">The lower limit of integration.</param>
        /// <param name="b">The upper limit of integration.</param>
        /// <param name="abserr">Estimate of the modulus of the absolute error, which should not exceed abs(i-result).</param>
        /// <param name="resabs">Approximation to the integral j.</param>
        /// <param name="resasc">Approximation to the integal of abs(f-i/(b-a)) over (a,b).</param>
        /// <returns>Approximation to the integral.</returns>
        public static double G_K15I(HardRealFunction f, double boun, int inf, double a, double b, ref double abserr,
                                    ref double resabs, ref double resasc)
        {
            double[] XGK15 = {
                                 0.99145537112081263921, 0.94910791234275852453, 0.86486442335976907279,
                                 0.74153118559939443986, 0.58608723546769113029, 0.40584515137739716691,
                                 0.20778495500789846760, 0.00000000000000000000
                             };
            double[] WGK15 = {
                                 0.02293532201052922496, 0.06309209262997855329, 0.10479001032225018384,
                                 0.14065325971552591875, 0.16900472663926790283, 0.19035057806478540991,
                                 0.20443294007529889241, 0.20948214108472782801
                             };
            double[] WG7 = {
                               0.12948496616886969327, 0.27970539148927666790, 0.38183005050511894495,
                               0.41795918367346938776
                           };

            double[] fv1 = new double[8];
            double[] fv2 = new double[8];
            double absc;
            double absc1;
            double absc2;
            double centr;
            double dinf;
            double fc;
            double fsum;
            double fval1;
            double fval2;
            double hlgth;
            double resg;
            double resk;
            double reskh;
            double result;
            double tabsc1;
            double tabsc2;
            int j;

            dinf = ((((1.0)) < (inf)) ? ((1.0)) : (inf));
            centr = 0.5 * (a + b);
            hlgth = 0.5 * (b - a);
            tabsc1 = boun + dinf * (1.0 - centr) / centr;
            fval1 = f.SolveAt(tabsc1);
            if (inf == 2)
            {
                fval1 += f.SolveAt(-tabsc1);
            }
            fc = (fval1 / centr) / centr;
            resg = fc * WG7[3];
            resk = fc * WGK15[7];
            resabs = Math.Abs(resk);
            for (j = 0; j < 7; j++)
            {
                absc = hlgth * XGK15[j];
                absc1 = centr - absc;
                absc2 = centr + absc;
                tabsc1 = boun + dinf * (1.0 - absc1) / absc1;
                tabsc2 = boun + dinf * (1.0 - absc2) / absc2;
                fval1 = f.SolveAt(tabsc1);
                fval2 = f.SolveAt(tabsc2);
                if (inf == 2)
                {
                    fval1 += f.SolveAt(-tabsc1);
                    fval2 += f.SolveAt(-tabsc2);
                }
                fval1 = (fval1 / absc1) / absc1;
                fval2 = (fval2 / absc2) / absc2;
                fv1[j] = fval1;
                fv2[j] = fval2;
                fsum = fval1 + fval2;
                if ((j & 1) != 0) // odd 'j's are truncated 
                {
                    resg += WG7[j / 2] * fsum;
                }
                resk += WGK15[j] * fsum;
                resabs = (resabs) + WGK15[j] * (Math.Abs(fval1) + Math.Abs(fval2));
            }
            reskh = resk * 0.5;
            resasc = WGK15[7] * Math.Abs(fc - reskh);
            for (j = 0; j < 7; j++)
            {
                resasc = (resasc) + WGK15[j] * (Math.Abs(fv1[j] - reskh) + Math.Abs(fv2[j] - reskh));
            }
            result = resk * hlgth;
            resabs = (resabs) * hlgth;
            resasc = (resasc) * hlgth;
            abserr = Math.Abs((resk - resg) * hlgth);
            if ((resasc != 0.0) && (abserr != 0.0))
            {
                abserr = (resasc) *
                         (((1.0) < (Math.Pow((200.0 * (abserr) / (resasc)), 1.5)))
                              ? (1.0)
                              : (Math.Pow((200.0 * (abserr) / (resasc)), 1.5)));
            }
            if (resabs > Double.MinValue / (50.0 * Double.Epsilon))
            {
                abserr = (((Double.Epsilon * 50.0 * (resabs)) > ((abserr)))
                              ? (Double.Epsilon * 50.0 * (resabs))
                              : ((abserr)));
            }
            return result;
        }
    }
}